module initInterpMindist

  ! ------------------------------------------------------------------------
  ! This module contains most of the "interesting" logic of initInterp, in terms of
  ! finding the input column (or landunit, patch, etc.) to use as a template for each
  ! output column (etc.).
  !
  ! This is in a separate module to facilitate unit testing, since the full initInterp
  ! involves some awkward dependencies.
  ! ------------------------------------------------------------------------

#include "shr_assert.h"
  use shr_kind_mod   , only: r8 => shr_kind_r8
  use shr_log_mod    , only : errMsg => shr_log_errMsg
  use clm_varctl     , only: iulog
  use abortutils     , only: endrun
  use spmdMod        , only: masterproc
  use clm_varcon     , only: spval, re
  use glcBehaviorMod , only: glc_behavior_type

  implicit none
  private
  save

  ! Public methods

  public :: set_mindist
  public :: set_single_match

  ! Public types

  type, public :: subgrid_special_indices_type
     integer :: ipft_not_vegetated
     integer :: icol_vegetated_or_bare_soil
     integer :: icol_urban_roof
     integer :: icol_urban_sunwall
     integer :: icol_urban_shadewall
     integer :: icol_urban_impervious_road
     integer :: icol_urban_pervious_road
     integer :: ilun_vegetated_or_bare_soil
     integer :: ilun_crop
     integer :: ilun_landice
     integer :: ilun_urban_TBD
     integer :: ilun_urban_HD
     integer :: ilun_urban_MD
   contains
     procedure :: is_vegetated_landunit  ! returns true if the given landunit type is natural veg or crop
     procedure :: is_urban_landunit      ! returns true if the given landunit type is urban
  end type subgrid_special_indices_type

  type, public :: subgrid_type
     character(len=16) :: name               ! pft, column, landunit, gridcell
     integer , pointer :: ptype(:) => null() ! used for patch type 
     integer , pointer :: ctype(:) => null() ! used for patch or col type
     integer , pointer :: ltype(:) => null() ! used for pft, col or lun type
     real(r8), pointer :: topoglc(:) => null()
     real(r8), pointer :: lat(:)
     real(r8), pointer :: lon(:)
     real(r8), pointer :: coslat(:)
   contains
     procedure :: print_point  ! print info about one point
     final :: destroy_subgrid_type
  end type subgrid_type

  ! Private methods

  private :: set_glc_must_be_same_type
  private :: set_ice_adjustable_type
  private :: do_fill_missing_with_natveg
  private :: do_fill_missing_urban_with_HD
  private :: is_sametype
  private :: is_baresoil
  private :: is_urban_HD

  character(len=*), parameter, private :: sourcefile = &
       __FILE__

contains

  !-----------------------------------------------------------------------
  subroutine print_point(this, index, unit)
    !
    ! !DESCRIPTION:
    ! Print info about one point in a subgrid_type object
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(subgrid_type), intent(in) :: this
    integer            , intent(in) :: index
    integer            , intent(in) :: unit ! unit to which we should write the info
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'print_point'
    !-----------------------------------------------------------------------

    write(unit,*) 'subgrid level, index = ',&
         this%name, index
    write(unit,*) 'lat, lon = ', this%lat(index), ', ', this%lon(index)
    if (associated(this%ltype)) then
       write(unit,*) 'ltype: ', this%ltype(index)
    end if
    if (associated(this%ctype)) then
       write(unit,*) 'ctype: ', this%ctype(index)
    end if
    if (associated(this%ptype)) then
       write(unit,*) 'ptype: ', this%ptype(index)
    end if
    
  end subroutine print_point

  !-----------------------------------------------------------------------
  subroutine destroy_subgrid_type(this)
    !
    ! !DESCRIPTION:
    ! Finalize routine for subgrid_type
    !
    ! !ARGUMENTS:
    type(subgrid_type) :: this
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'destroy_subgrid_type'
    !-----------------------------------------------------------------------

    if (associated(this%ptype)) then
       deallocate(this%ptype)
    end if

    if (associated(this%ctype)) then
       deallocate(this%ctype)
    end if

    if (associated(this%ltype)) then
       deallocate(this%ltype)
    end if

    if (associated(this%topoglc)) then
       deallocate(this%topoglc)
    end if

    if (associated(this%lat)) then
       deallocate(this%lat)
    end if

    if (associated(this%lon)) then
       deallocate(this%lon)
    end if

    if (associated(this%coslat)) then
       deallocate(this%coslat)
    end if

  end subroutine destroy_subgrid_type

  !=======================================================================

  subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgrido, &
       subgrid_special_indices, glc_behavior, glc_elevclasses_same, &
       fill_missing_with_natveg, fill_missing_urban_with_HD, mindist_index)

    ! --------------------------------------------------------------------
    ! arguments
    integer            , intent(in)  :: begi, endi 
    integer            , intent(in)  :: bego, endo 
    logical            , intent(in)  :: activei(begi:endi) 
    logical            , intent(in)  :: activeo(bego:endo) 
    type(subgrid_type) , intent(in)  :: subgridi
    type(subgrid_type) , intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    type(glc_behavior_type), intent(in) :: glc_behavior

    ! True if the number and bounds of glacier elevation classes are the same in the
    ! input and output; false otherwise.
    logical            , intent(in)  :: glc_elevclasses_same

    ! If false: if an output type cannot be found in the input, code aborts
    ! If true: if a non-urban output type cannot be found in the input, fill with closest natural
    ! veg column (using bare soil for patch-level variables)
    !
    ! NOTE: always treated as true for natural veg and crop landunits/columns/patches in
    ! the output - e.g., if we can't find the right column type to fill crop, we always
    ! use the closest natural veg column, regardless of the value of this flag.
    logical            , intent(in)  :: fill_missing_with_natveg


    ! If false: if an urban output type cannot be found in the input, code aborts
    ! If true: if an urban output type cannot be found in the input, fill with closest urban HD
    logical            , intent(in)  :: fill_missing_urban_with_HD

    integer            , intent(out) :: mindist_index(bego:endo)
    !
    ! local variables
    real(r8) :: dx,dy
    real(r8) :: distmin,dist,hgtdiffmin,hgtdiff    
    integer  :: nsizei, nsizeo
    integer  :: ni,no,nmin,ier,n,noloc
    logical  :: topoglc_present
    logical  :: closest

    ! Whether two glc columns/patches must be the same column/patch type to be
    ! considered the same type. This is only valid for glc points, and is only valid
    ! for subgrid name = 'pft' or 'column'.
    logical  :: glc_must_be_same_type_o(bego:endo)

    character(len=*), parameter :: subname = 'set_mindist'
    ! --------------------------------------------------------------------

    if (associated(subgridi%topoglc) .and. associated(subgrido%topoglc)) then
       topoglc_present = .true.
    else
       topoglc_present = .false.
    end if

    mindist_index(bego:endo) = 0
    distmin = spval

    call set_glc_must_be_same_type(bego=bego, endo=endo, dimname=subgrido%name, &
         glc_elevclasses_same = glc_elevclasses_same, glc_behavior=glc_behavior, &
         glc_must_be_same_type_o=glc_must_be_same_type_o)

!$OMP PARALLEL DO PRIVATE (ni,no,n,nmin,distmin,dx,dy,dist,closest,hgtdiffmin,hgtdiff)
    do no = bego,endo

       ! Only interpolate onto active points. Otherwise, the mere act of running
       ! init_interp (e.g., of a file onto itself) would lead to changes in a bunch of
       ! inactive points - e.g., going from their cold start initial conditions to some
       ! spunup initial conditions (from the closest active point of that type). This
       ! could potentially lead to different behavior in a transient run, if those points
       ! later became active; that's undesirable.
       if (activeo(no)) then 
          nmin    = 0
          distmin = spval
          hgtdiffmin = spval
          do ni = begi,endi
             if (activei(ni)) then
                if (is_sametype(ni = ni, no = no, &
                     subgridi = subgridi, subgrido = subgrido, &
                     subgrid_special_indices = subgrid_special_indices, &
                     glc_must_be_same_type = glc_must_be_same_type_o(no), &
                     veg_patch_just_considers_ptype = .true., &
                     do_fill_missing_urban_with_HD = .false.)) then
                   dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re
                   dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * &
                        0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni))
                   dist = dx*dx + dy*dy
                   if (topoglc_present) then
                      hgtdiff = abs(subgridi%topoglc(ni) - subgrido%topoglc(no))
                   end if
                   closest = .false.
                   if ( dist < distmin ) then
                      closest = .true.
                      distmin = dist
                      nmin = ni
                      if (topoglc_present) then
                         hgtdiffmin = hgtdiff
                      end if
                   end if
                   if (.not. closest) then
                      ! In *some* cases: For glacier points, we first find the closest
                      ! point in lat-lon space, without consideration for column or patch
                      ! type (above). Then, within that closest point, we find the closest
                      ! column in topographic space; this second piece is done here. Note
                      ! that this ordering means that we could choose a column with a very
                      ! different topographic height from the target column, if it is
                      ! closer in lat-lon space than any similar-height columns.
                      if (topoglc_present) then
                         hgtdiff = abs(subgridi%topoglc(ni) - subgrido%topoglc(no))
                         if ((dist == distmin) .and. (hgtdiff < hgtdiffmin)) then
                            closest = .true.
                            hgtdiffmin = hgtdiff
                            distmin = dist
                            nmin = ni
                         end if
                      end if
                   end if
                end if
             end if
          end do
          
          ! Note that do_fill_missing_with_natveg below will return .false. for pfts and columnns associated
          ! with urban landunits so that the fill missing with bare soil will be implemented only for
          ! non-urban types (pfts, columns, landunits, gridcells).

          ! If non-urban output type is not contained in input dataset, then use closest bare soil,
          ! if this point is one for which we fill missing with natveg.
          if ( distmin == spval .and. &
               do_fill_missing_with_natveg( &
               fill_missing_with_natveg, no, subgrido, subgrid_special_indices)) then
             do ni = begi, endi
                if (activei(ni)) then
                   if ( is_baresoil(ni, subgridi, subgrid_special_indices)) then
                      dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re
                      dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * &
                           0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni))
                      dist = dx*dx + dy*dy
                      if ( dist < distmin )then
                         distmin = dist
                         nmin = ni
                      end if
                   end if
                end if
             end do

          ! If urban output type is not contained in input dataset, then use closest urban HD,
          ! if this point is one for which we fill missing urban with urban HD.
          else if (distmin == spval &
               .and. do_fill_missing_urban_with_HD( &
               fill_missing_urban_with_HD, no, subgrido, subgrid_special_indices)) then
             do ni = begi, endi
                if (activei(ni)) then
                   ! We need to call is_sametype for pfts and columns here to make sure that each
                   ! urban input pft and column type matches the output pft and column type. We don't
                   ! want to call it for landunits because they intentionally won't be the same type
                   ! (since we are filling missing urban landunits with HD)
                   if (subgrido%name .eq. 'landunit') then
                      if ( is_urban_HD(ni, subgridi, subgrid_special_indices)) then
                         dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re
                         dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * &
                              0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni))
                         dist = dx*dx + dy*dy
                         if ( dist < distmin )then
                            distmin = dist
                            nmin = ni
                         end if
                      end if
                   else
                      if (is_sametype(ni = ni, no = no, &
                          subgridi = subgridi, subgrido = subgrido, &
                          subgrid_special_indices = subgrid_special_indices, &
                          glc_must_be_same_type = glc_must_be_same_type_o(no), &
                          veg_patch_just_considers_ptype = .false., &
                          do_fill_missing_urban_with_HD = .true.)) then
                         if ( is_urban_HD(ni, subgridi, subgrid_special_indices)) then
                            dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re
                            dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * &
                                 0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni))
                            dist = dx*dx + dy*dy
                            if ( dist < distmin )then
                               distmin = dist
                               nmin = ni
                            end if
                         end if
                      end if
                   end if
                end if
             end do
          end if

          ! Error conditions
          if ( distmin == spval )then
             write(iulog,*) 'ERROR initInterp set_mindist: &
                  &Cannot find any input points matching output point:'
             call subgrido%print_point(no, iulog)
             write(iulog,*) ' '
             write(iulog,*) 'If this is an urban type'
             write(iulog,*) '(ltype = ', subgrid_special_indices%ilun_urban_TBD, &
                            ',', subgrid_special_indices%ilun_urban_HD, &
                            ', or', subgrid_special_indices%ilun_urban_MD, ')'
             write(iulog,*) 'then consider rerunning with the following in user_nl_clm:'
             write(iulog,*) 'init_interp_fill_missing_urban_with_HD = .true.'
             write(iulog,*) 'However, note that this will fill all urban missing types in the output'
             write(iulog,*) 'with the closest urban high density (HD) type in the input'
             write(iulog,*) 'So, you should consider whether that is what you want.'
             write(iulog,*) ' '
             write(iulog,*) 'If this is a non-urban type'
             write(iulog,*) '(ltype \= ',subgrid_special_indices%ilun_urban_TBD, &
                            ',', subgrid_special_indices%ilun_urban_HD, &
                            ', or', subgrid_special_indices%ilun_urban_MD, ')'
             write(iulog,*) 'consider rerunning with the following in user_nl_clm:'
             write(iulog,*) 'init_interp_fill_missing_with_natveg = .true.'
             write(iulog,*) 'However, note that this will fill all non-urban missing types in the output'
             write(iulog,*) 'with the closest natural veg column in the input'
             write(iulog,*) '(using bare soil for patch-level variables).'
             write(iulog,*) 'So, you should consider whether that is what you want.'
             write(iulog,*) errMsg(sourcefile, __LINE__)
             call endrun(msg=subname// &
                  ' ERROR: Cannot find any input points matching output point')
          end if

          mindist_index(no) = nmin

       end if ! end if activeo block
    end do
!$OMP END PARALLEL DO
    
  end subroutine set_mindist

  !-----------------------------------------------------------------------
  subroutine set_single_match(begi, endi, bego, endo, activeo, subgridi, subgrido, &
       subgrid_special_indices, glc_behavior, glc_elevclasses_same, &
       mindist_index)
    !
    ! !DESCRIPTION:
    ! This is an alternative to set_mindist that can work when there is a single matching
    ! input point for each output point (or possibly 0 matching inputs for an inactive
    ! output point). In this method, we only look for points that share the same lat/lon
    ! as the output point.
    !
    ! This method is intended to be used for the case where we want to copy areas from the
    ! input file to the output file (rather than maintaining areas from the surface
    ! dataset); some of the logic here is set up the way it is in order to support this
    ! copying of areas.
    !
    ! !ARGUMENTS:
    integer            , intent(in)  :: begi, endi
    integer            , intent(in)  :: bego, endo
    logical            , intent(in)  :: activeo(bego:endo)
    type(subgrid_type) , intent(in)  :: subgridi
    type(subgrid_type) , intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    type(glc_behavior_type), intent(in) :: glc_behavior

    ! True if the number and bounds of glacier elevation classes are the same in the
    ! input and output; false otherwise.
    !
    ! Must be true for this method! (Otherwise we may not be able to find a single unique
    ! input point for each output point.)
    logical            , intent(in)  :: glc_elevclasses_same

    integer            , intent(out) :: mindist_index(bego:endo)
    !
    ! !LOCAL VARIABLES:
    integer  :: ni, no, ni_match
    logical  :: found
    real(r8) :: dx, dy
    logical  :: ni_sametype

    ! Whether two glc columns/patches must be the same column/patch type to be
    ! considered the same type. This is only valid for glc points, and is only valid
    ! for subgrid name = 'pft' or 'column'.
    logical  :: glc_must_be_same_type_o(bego:endo)

    ! Tolerance in lat/lon for considering a point to be at the same location
    real(r8) :: same_point_tol = 1.e-14_r8

    character(len=*), parameter :: subname = 'set_single_match'
    !-----------------------------------------------------------------------

    if (.not. glc_elevclasses_same) then
       write(iulog,*) subname// &
            ' ERROR: For this init_interp method, the number and bounds of'
       write(iulog,*) 'glacier elevation classes must be the same in input and output'
       call endrun(msg=subname//' ERROR: glc_elevclasses_same must be true for this method')
    end if

    call set_glc_must_be_same_type(bego=bego, endo=endo, dimname=subgrido%name, &
         glc_elevclasses_same = glc_elevclasses_same, glc_behavior=glc_behavior, &
         glc_must_be_same_type_o=glc_must_be_same_type_o)

!$OMP PARALLEL DO PRIVATE (ni,no,ni_match,found,dx,dy,ni_sametype)
    do no = bego, endo
       ni_match = 0
       found = .false.
       do ni = begi, endi
          dx = abs(subgrido%lon(no)-subgridi%lon(ni))
          dy = abs(subgrido%lat(no)-subgridi%lat(ni))
          if (dx < same_point_tol .and. dy < same_point_tol) then
             ni_sametype = is_sametype(ni = ni, no = no, &
                  subgridi = subgridi, subgrido = subgrido, &
                  subgrid_special_indices = subgrid_special_indices, &
                  glc_must_be_same_type = glc_must_be_same_type_o(no), &
                  veg_patch_just_considers_ptype = .false., &
                  do_fill_missing_urban_with_HD = .false.)
             if (ni_sametype) then
                if (found) then
                   write(iulog,*) subname// &
                        ' ERROR: found multiple input points matching output point'
                   call subgrido%print_point(no, iulog)
                   write(iulog,*) ' '
                   write(iulog,*) 'For this init_interp_method: for a given output point,'
                   write(iulog,*) 'we expect to find exactly one input point at the same'
                   write(iulog,*) 'location with the same type.'
                   write(iulog,*) '(This most likely indicates a problem in CTSM, such as'
                   write(iulog,*) 'having multiple columns with the same column type on a'
                   write(iulog,*) 'given gridcell.)'
                   call endrun(msg=subname// &
                        ' ERROR: found multiple input points matching output point')
                else
                   found = .true.
                   ni_match = ni
                end if
             end if
          end if
       end do

       if (found) then
          mindist_index(no) = ni_match
       else
          mindist_index(no) = 0

          if (activeo(no)) then
             write(iulog,*) subname// &
                  ' ERROR: cannot find any input points matching output point'
             call subgrido%print_point(no, iulog)
             write(iulog,*) ' '
             write(iulog,*) 'For this init_interp_method: for a given active output point,'
             write(iulog,*) 'we expect to find exactly one input point at the same'
             write(iulog,*) 'location with the same type.'
             write(iulog,*) 'Note that this requires the input and output grids to be identical.'
             write(iulog,*) '(If you need to interpolate to a different resolution, then'
             write(iulog,*) 'you need to use a different init_interp_method.)'
             call endrun(msg=subname//' ERROR: cannot find any input points matching output point')
          end if
       end if
    end do
!$OMP END PARALLEL DO

  end subroutine set_single_match

  !-----------------------------------------------------------------------
  subroutine set_glc_must_be_same_type(bego, endo, dimname, &
       glc_elevclasses_same, glc_behavior, &
       glc_must_be_same_type_o)
    !
    ! !DESCRIPTION:
    ! Sets the glc_must_be_same_type_o array for each output ice point
    !
    ! This array will be set to true for ice output columns/patches for which the
    ! column/patch type must match the input column/patch type to be considered the same
    ! type. This is only valid for ice points, and is only valid for dimname = 'pft' or
    ! 'column' - for others, the value is undefined.
    !
    ! This assumes that bego and endo match the bounds that are used elsewhere in the
    ! model - i.e., for the decomposition within init_interp to match the model's
    ! decomposition.
    !
    ! !ARGUMENTS:
    integer                 , intent(in)  :: bego    ! beginning index for output points
    integer                 , intent(in)  :: endo    ! ending index for output points
    character(len=*)        , intent(in)  :: dimname ! 'pft', 'column', etc.

    ! True if the number and bounds of glacier elevation classes are the same in the
    ! input and output; false otherwise.
    logical, intent(in) :: glc_elevclasses_same

    type(glc_behavior_type), intent(in)  :: glc_behavior

    logical, intent(out)  :: glc_must_be_same_type_o( bego: ) ! see description above

    !
    ! !LOCAL VARIABLES:
    integer :: no

    ! Whether each output ice point has adjustable column type; only valid for ice
    ! points, and only valid for subgrid name = pft or column
    logical  :: ice_adjustable_type_o(bego:endo)

    character(len=*), parameter :: subname = 'set_glc_must_be_same_type'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(glc_must_be_same_type_o) == (/endo/)), sourcefile, __LINE__)

    if (.not. glc_elevclasses_same) then
       ! If the number or bounds of the elevation classes differ between input and
       ! output, then we ignore the column and patch types of glc points when
       ! looking for the same type - instead using closest topographic height as a
       ! tie-breaker between equidistant columns/patches.
       glc_must_be_same_type_o(bego:endo) = .false.
    else
       call set_ice_adjustable_type(bego=bego, endo=endo, dimname=dimname, &
            glc_behavior=glc_behavior, ice_adjustable_type_o=ice_adjustable_type_o)
       do no = bego, endo
          if (ice_adjustable_type_o(no)) then
             ! If glc points in this output cell have adjustable type, then we ignore
             ! the column and patch types of glc points when looking for the same
             ! type: we want to find the closest glc point in lat-lon space without
             ! regards for column/patch type, because the column/patch types may change
             ! at runtime.
             glc_must_be_same_type_o(no) = .false.
          else
             ! Otherwise, we require the column and patch types to be the same between
             ! input and output for this glc output point, as is the case for most
             ! other landunits. This is important for a case with interpolation to give
             ! bit-for-bit answers with a case without interpolation (since glc
             ! topographic heights can change after initialization, so we can't always
             ! rely on the point with closest topographic height to be the "right" one to
             ! pick as the source for interpolation).
             glc_must_be_same_type_o(no) = .true.
          end if
       end do
    end if

  end subroutine set_glc_must_be_same_type

  !-----------------------------------------------------------------------
  subroutine set_ice_adjustable_type(bego, endo, dimname, glc_behavior, &
       ice_adjustable_type_o)
    !
    ! !DESCRIPTION:
    ! Sets the ice_adjustable_type_o array for each output ice point
    !
    ! This array will be set to true for ice points that have adjustable column type
    ! and false for ice points that do not. The value will be undefined for non-ice
    ! points.
    !
    ! This can only be called for the output, not the input!
    !
    ! The output array is only valid for dimname = pft or column; for others, the output
    ! values are undefined.
    !
    ! This assumes that bego and endo match the bounds that are used elsewhere in the
    ! model - i.e., for the decomposition within init_interp to match the model's
    ! decomposition.
    !
    ! !ARGUMENTS:
    integer                 , intent(in)  :: bego    ! beginning index for output points
    integer                 , intent(in)  :: endo    ! ending index for output points
    character(len=*)        , intent(in)  :: dimname ! 'pft', 'column', etc.
    type(glc_behavior_type) , intent(in)  :: glc_behavior
    logical                 , intent(out) :: ice_adjustable_type_o( bego: ) ! see documentation above
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'set_ice_adjustable_type'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(ice_adjustable_type_o) == (/endo/)), sourcefile, __LINE__)

    select case (dimname)
    case ('pft')
       call glc_behavior%patches_have_dynamic_type_array(bego, endo, &
            ice_adjustable_type_o(bego:endo))
    case ('column')
       call glc_behavior%cols_have_dynamic_type_array(bego, endo, &
            ice_adjustable_type_o(bego:endo))
    case ('landunit', 'gridcell')
       ! Do nothing: ice_adjustable_type_o will be left undefined
    case default
       call endrun(subname//' ERROR: unexpected dimname: '//trim(dimname)//&
            errMsg(sourcefile, __LINE__))
    end select

  end subroutine set_ice_adjustable_type

  !-----------------------------------------------------------------------
  function do_fill_missing_with_natveg(fill_missing_with_natveg, &
       no, subgrido, subgrid_special_indices)
    !
    ! !DESCRIPTION:
    ! Returns true if the given non-urban output point, if missing, should be filled with the
    ! closest natural veg point.
    !
    ! !ARGUMENTS:
    logical :: do_fill_missing_with_natveg  ! function result

    ! whether we should fill ALL missing points with natveg
    logical, intent(in) :: fill_missing_with_natveg

    integer           , intent(in)  :: no
    type(subgrid_type), intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'do_fill_missing_with_natveg'
    !-----------------------------------------------------------------------

    if (subgrido%name == 'gridcell') then
       ! It makes no sense to try to fill missing with natveg for gridcell-level values
       do_fill_missing_with_natveg = .false.
    else if (fill_missing_with_natveg .and. .not. subgrid_special_indices%is_urban_landunit(subgrido%ltype(no))) then
       ! User has asked for all non-urban missing points to be filled with natveg
       do_fill_missing_with_natveg = .true.
    else if (subgrid_special_indices%is_vegetated_landunit(subgrido%ltype(no))) then
       ! Even if user hasn't asked for it, we fill missing vegetated points (natural veg
       ! and crop) with the closest natveg point. This is mainly to support the common
       ! use case of interpolating non-crop to crop, but also supports adding a new PFT
       ! type.
       do_fill_missing_with_natveg = .true.
    else
       do_fill_missing_with_natveg = .false.
    end if

  end function do_fill_missing_with_natveg

  !-----------------------------------------------------------------------
  function do_fill_missing_urban_with_HD(fill_missing_urban_with_HD, &
       no, subgrido, subgrid_special_indices)
    !
    ! !DESCRIPTION:
    ! Returns true if the given urban output point, if missing, should be filled with the
    ! closest urban HD point.
    !
    ! !ARGUMENTS:
    logical :: do_fill_missing_urban_with_HD  ! function result

    ! whether we should fill ALL missing points with urban HD
    logical, intent(in) :: fill_missing_urban_with_HD

    integer           , intent(in)  :: no
    type(subgrid_type), intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'do_fill_missing_urban_with_HD'
    !-----------------------------------------------------------------------

    if (subgrido%name == 'gridcell') then
       ! It makes no sense to try to fill missing with urban HD for gridcell-level values
       do_fill_missing_urban_with_HD = .false.
    else if (fill_missing_urban_with_HD) then
       ! User has asked for all missing urban points to be filled with urban HD
       do_fill_missing_urban_with_HD = .true.
    else
       do_fill_missing_urban_with_HD = .false.
    end if

  end function do_fill_missing_urban_with_HD

  !=======================================================================

  logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indices, &
                  glc_must_be_same_type, veg_patch_just_considers_ptype, &
                  do_fill_missing_urban_with_HD)

    ! --------------------------------------------------------------------
    ! arguments
    integer           , intent(in)  :: ni 
    integer           , intent(in)  :: no 
    type(subgrid_type), intent(in)  :: subgridi
    type(subgrid_type), intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices

    ! Whether two glc columns/patches must be the same column/patch type to be
    ! considered the same type. This is only valid for glc points, and is only valid
    ! for subgrid name = 'pft' or 'column'.
    logical, intent(in) :: glc_must_be_same_type

    ! For vegetated patches (natural veg or crop): if veg_patch_just_considers_ptype is
    ! true, then we consider two vegetated patches to be the same type if they have the
    ! same ptype, without regard to their column and landunit types (as long as both the
    ! input and the output are on either the natural veg or crop landunit). This is
    ! needed to handle the generic crop properly when interpolating from non-crop to
    ! crop, or vice versa.
    !
    ! If false, then they need to have the same column and landunit types, too (as is the
    ! general case).
    logical, intent(in) :: veg_patch_just_considers_ptype

    ! If True, we allow for landunits to be different when checking if pft and column are
    ! the same type, to allow for HD fill of missing urban output points.
    logical, intent(in) :: do_fill_missing_urban_with_HD

    ! For urban columns/patches
    ! --------------------------------------------------------------------

    is_sametype = .false.

    if (trim(subgridi%name) == 'pft' .and. trim(subgrido%name) == 'pft') then
       if ( .not. glc_must_be_same_type .and. &
            subgridi%ltype(ni) == subgrid_special_indices%ilun_landice .and. &
            subgrido%ltype(no) == subgrid_special_indices%ilun_landice) then
          is_sametype = .true.
       else if (veg_patch_just_considers_ptype .and. &
            subgrid_special_indices%is_vegetated_landunit(subgrido%ltype(no))) then
          ! See comment attached to the declaration of veg_patch_just_considers_ptype for
          ! rationale for this logic.
          !
          ! TODO(wjs, 2015-09-15) If we ever allow the same PFT to appear on multiple
          ! columns within a given grid cell, then this logic will need to be made
          ! somewhat more complex: e.g., preferably take something from the same column
          ! type, but if we can't find anything from the same column type, then ignore
          ! column type.

          if (subgrid_special_indices%is_vegetated_landunit(subgridi%ltype(ni)) .and. &
               subgridi%ptype(ni) == subgrido%ptype(no)) then
             is_sametype = .true.
          end if
       else if (subgridi%ptype(ni) == subgrido%ptype(no) .and. &
                subgridi%ctype(ni) == subgrido%ctype(no) .and. &
                do_fill_missing_urban_with_HD) then
          is_sametype = .true.
       else if (subgridi%ptype(ni) == subgrido%ptype(no) .and. &
                subgridi%ctype(ni) == subgrido%ctype(no) .and. &
                subgridi%ltype(ni) == subgrido%ltype(no)) then
          is_sametype = .true.
       end if
    else if (trim(subgridi%name) == 'column' .and. trim(subgrido%name) == 'column') then
       if ( .not. glc_must_be_same_type .and. &
            subgridi%ltype(ni) == subgrid_special_indices%ilun_landice  .and. &
            subgrido%ltype(no) == subgrid_special_indices%ilun_landice ) then
          is_sametype = .true.
       else if (subgridi%ctype(ni) == subgrido%ctype(no) .and. &
                do_fill_missing_urban_with_HD) then
          is_sametype = .true.
       else if (subgridi%ctype(ni) == subgrido%ctype(no) .and. &
                subgridi%ltype(ni) == subgrido%ltype(no)) then
          is_sametype = .true.
       end if
    else if (trim(subgridi%name) == 'landunit' .and. trim(subgrido%name) == 'landunit') then
       if (subgridi%ltype(ni) == subgrido%ltype(no)) then
          is_sametype = .true.
       end if
    else if (trim(subgridi%name) == 'gridcell' .and. trim(subgrido%name) == 'gridcell') then
       is_sametype = .true.
    else 
       if (masterproc) then
          write(iulog,*)'ERROR interpinic: is_sametype check on input and output type not supported'
          write(iulog,*)'typei = ',trim(subgridi%name)
          write(iulog,*)'typeo = ',trim(subgrido%name)
       end if
       call endrun(msg=errMsg(sourcefile, __LINE__))
    end if

  end function is_sametype

  !=======================================================================

  logical function is_baresoil (n, subgrid, subgrid_special_indices)

    ! --------------------------------------------------------------------
    ! arguments
    integer           , intent(in)  :: n 
    type(subgrid_type), intent(in)  :: subgrid
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    ! --------------------------------------------------------------------

    is_baresoil = .false.

    if (subgrid%name == 'pft') then
       if (subgrid%ptype(n) == subgrid_special_indices%ipft_not_vegetated .and. &
            subgrid%ctype(n) == subgrid_special_indices%icol_vegetated_or_bare_soil .and. &
            subgrid%ltype(n) == subgrid_special_indices%ilun_vegetated_or_bare_soil) then
          is_baresoil = .true.
       end if
    else if (subgrid%name == 'column') then
       if (subgrid%ctype(n) == subgrid_special_indices%icol_vegetated_or_bare_soil .and. &
            subgrid%ltype(n) == subgrid_special_indices%ilun_vegetated_or_bare_soil) then
          is_baresoil = .true.
       end if
    else if (subgrid%name == 'landunit') then
       if (subgrid%ltype(n) == subgrid_special_indices%ilun_vegetated_or_bare_soil) then
          is_baresoil = .true.
       end if
    else 
       if (masterproc) then
          write(iulog,*)'ERROR interpinic: is_baresoil subgrid type ',subgrid%name,' not supported'
       end if
       call endrun(msg=errMsg(sourcefile, __LINE__))
    end if

  end function is_baresoil

  !-----------------------------------------------------------------------
  logical function is_urban_HD (n, subgrid, subgrid_special_indices)

    ! --------------------------------------------------------------------
    ! arguments
    integer           , intent(in)  :: n
    type(subgrid_type), intent(in)  :: subgrid
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    ! --------------------------------------------------------------------

    is_urban_HD = .false.

    if (subgrid%name == 'pft' .or. subgrid%name == 'column' .or. subgrid%name == 'landunit') then
       if (subgrid%ltype(n) == subgrid_special_indices%ilun_urban_HD) then
          is_urban_HD = .true.
       end if
    else
       if (masterproc) then
          write(iulog,*)'ERROR interpinic: is_urban_HD subgrid type ',subgrid%name,' not supported'
       end if
       call endrun(msg=errMsg(sourcefile, __LINE__))
    end if

  end function is_urban_HD

  !-----------------------------------------------------------------------
  function is_vegetated_landunit(this, ltype)
    !
    ! !DESCRIPTION:
    ! Returns true if the given landunit type is vegetated: either natural veg or crop
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    logical :: is_vegetated_landunit  ! function result
    class(subgrid_special_indices_type), intent(in) :: this
    integer, intent(in) :: ltype   ! landunit type of interest
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'is_vegetated_landunit'
    !-----------------------------------------------------------------------

    if (ltype == this%ilun_vegetated_or_bare_soil .or. &
         ltype == this%ilun_crop) then
       is_vegetated_landunit = .true.
    else
       is_vegetated_landunit = .false.
    end if

  end function is_vegetated_landunit

  function is_urban_landunit(this, ltype)
    !
    ! !DESCRIPTION:
    ! Returns true if the given landunit type is urban
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    logical :: is_urban_landunit  ! function result
    class(subgrid_special_indices_type), intent(in) :: this
    integer, intent(in) :: ltype   ! landunit type of interest
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'is_urban_landunit'
    !-----------------------------------------------------------------------

    if (ltype == this%ilun_urban_TBD .or. ltype == this%ilun_urban_HD &
            .or. ltype == this%ilun_urban_MD) then
       is_urban_landunit = .true.
    else
       is_urban_landunit = .false.
    end if

  end function is_urban_landunit

end module initInterpMindist
